Exam 1 - Module 3
CQF Jan 2023

Amir Azmi (amirmazmi@gmail.com)


Question 1¶

Given the follow set of sample data, compare the different schemes for Geometric Brownian Motion (GBM)

$$ \begin{align} \text{Today’s stock price, } S_0 &= 100 \\[5pt] \text{Time, } T &= 1 \text{ year} \\[5pt] \text{volatility, } \sigma &= 20\% \\[5pt] \text{constant risk-free interest rate, } r &= 5\% \\[5pt] \end{align} $$

1.1 Different schemes for GBM¶

The following schemes will be considered for comparison

  • Euler-Maruyama scheme
  • Milstein scheme
  • Closed form solution

Begin by looking at the closed form solution.

1.1.1 Closed form solution¶

$$ S_{t + \delta t} = S_t \; exp \left\{ \; ( r - \frac{1}{2} \sigma^2 ) \; \delta t \;+\; \sigma \phi \sqrt{\delta t} \; \right\} \qquad \text{where} \qquad \phi \sim N(0,1) $$
Dr Ahmad, R. (2023). Module 1 Lecture 5 - "Simulating and Manipulating Stochastic Differential Equations"[PowerPoint slides].
CQF Program January 2023

1.1.2 Euler-Maruyama¶

The Euler-Maruyama scheme is an approximation by taking a Taylor Series expansion on the exact solution

Let,
$$ x = ( r - \frac{1}{2} \sigma^2 ) \; \delta t \;+\; \sigma \phi \sqrt{\delta t} \\ $$

Then Taylor Series expansion on the exponential with $a=0$,

$$ \begin{align} e^x &\approx e^a \left( 1 + (x-a) + \frac{1}{2} (x-a)^2 + \frac{1}{6} (x-a)^3 \;+\; ... \right) \\ &\approx 1 + x + \frac{1}{2} x^2 + \frac{1}{6} x^3 \;+\; ... \\ \end{align}\\ $$
Therefore, $$ exp \left\{ \; ( r - \frac{1}{2} \sigma^2 ) \; \delta t \;+\; \sigma \phi \sqrt{\delta t} \right\} \; \sim \; 1 \;+\; ( r - \frac{1}{2} \sigma^2 ) \; \delta t \;+\; \sigma \phi \sqrt{\delta t} \;+\; O(\delta t) \\ $$

Then dropping the higher order terms $O(\delta t)$,

$$ S_{t + \delta t} \;\sim\; S_t \; ( 1 \;+\; ( \; r - \frac{1}{2} \sigma^2 ) \; \delta t \;+\; \sigma \phi \sqrt{\delta t} \; ) \\ $$

Results in an approximation with an accuracy of the order $O(\sqrt{\delta t})$


Dr Ahmad, R. (2023). Module 1 Lecture 5 - "Simulating and Manipulating Stochastic Differential Equations" [PowerPoint slides].
CQF Program January 2023



1.1.3 Milstein Scheme¶

The Milstein Scheme, is an inclusion of one higher order term from the Taylor Expansion done in the Euler-Maruyama scheme.

$$ exp \left\{ \; ( r - \frac{1}{2} \sigma^2 ) \; \delta t \;+\; \sigma \phi \sqrt{\delta t} \right\} \; \sim \; 1 \;+\; ( r - \frac{1}{2} \sigma^2 ) \; \delta t \;+\; \sigma \phi \sqrt{\delta t} + \underbrace{ \frac{1}{2} \sigma^2 \phi^2 \delta t }_{ \text{Milstein correction}} + O(\delta t^2) \\ $$

Similarly, dropping the higher order terms $O(\delta t^2)$,

$$ S_{t + \delta t} \;\sim\; S_t \; \left( 1 \;+\; (\; r - \frac{1}{2} \sigma^2 ) \; \delta t \;+\; \sigma \phi \sqrt{\delta t} \;+\; \frac{1}{2} \sigma^2 \phi^2 \delta t \; \right) \\ \;\sim\; S_t \; \left( 1 \;+\; r \; \delta t \;+\; \sigma \phi \sqrt{\delta t} \;+\; \frac{1}{2} \sigma^2 ( \phi^2 -1) \; \delta t \; \right) \\ $$

Results in an approximation with an improved accuracy of the order $O(\delta t)$.


Dr Ahmad, R. (2023). Monte Carlo Tutorial - "Monte Carlo" [PowerPoint slides].
CQF Program January 2023

1.1.4 Monte Carlo simulation¶

To run the simulation, make the following assumptions:

  1. 1 year is 252 trading days, $\sigma$ is scaled by $\frac{1}{\sqrt{252}}$
  2. 1 day timesteps which means there are 252 timesteps ( $0 < t \leq T$ )
In [1]:
import numpy as np
import pandas as pd
import sympy as sp
sp.init_printing()
from functools import partial
from copy import deepcopy
from time import perf_counter
from IPython.display import display_html 
import plotly.express as px
import plotly.graph_objects as go
import plotly.io as pio
pio.renderers.default = 'notebook' #_connected'
from plotly.offline import init_notebook_mode
init_notebook_mode(connected=True)

from watermark import watermark
%load_ext watermark
%watermark -u -d -v -m -iv
Last updated: 2023-04-27

Python implementation: CPython
Python version       : 3.9.16
IPython version      : 8.12.0

Compiler    : GCC 7.5.0
OS          : Linux
Release     : 5.4.0-147-generic
Machine     : x86_64
Processor   : x86_64
CPU cores   : 8
Architecture: 64bit

sympy : 1.11.1
numpy : 1.24.2
pandas: 2.0.0
plotly: 5.14.1

Define the known constants

In [2]:
# define constants
S0 = 100
rate = 0.05
sigma = 0.2
T = 1
t = 252
delta_t = T/t
runs = 100_000

Initialize arrays for stock price and values for the random variables for repeatability across different schemes.

In [3]:
# initialize empty stock price array
S = np.zeros((t, runs))
S[0] = S0

# initialize random values
np.random.seed(20230419)
arr_rand = np.random.standard_normal(t*runs).reshape(t,runs)

Create a function for each scheme.

In [4]:
def exact(s_prev: float, r: float, vol: float, dt: float, randvar: int) -> np.array:
    return s_prev * np.exp( (r - 0.5 * vol**2) * dt + vol * np.sqrt(dt) * randvar)

def euler_maruyama(s_prev: float, r: float, vol: float, dt: float, randvar: int) -> np.array:
    return s_prev * ( 1 + r * dt + vol * np.sqrt(dt) * randvar)
    
def milstein(s_prev: float, r: float, vol: float, dt: float, randvar: int) -> np.array:
    return s_prev * ( 1 + r * dt + vol * np.sqrt(dt) * randvar + 0.5 * vol**2 * (randvar**2 - 1) * np.sqrt(dt))

ls_scheme = {'exact': partial(exact, r=rate, vol=sigma),
             'euler_maruyama': partial(euler_maruyama, r=rate, vol=sigma),
             'milstein': partial(milstein, r=rate, vol=sigma) }

Create a wrapper function to run Monte Carlo simulation for each scheme using the same random array.

In [5]:
def simulate_path(approx_func, dt: float, arr_price: np.array, arr_random: np.array) -> np.array:
    arr_out = deepcopy(arr_price)
    for i in range(0, arr_out.shape[0]-1):
        arr_out[i+1] = approx_func( s_prev=arr_out[i], dt=dt, randvar=arr_random[i] )
    return arr_out  
In [6]:
# helper function for long table display
def style_table(df: pd.DataFrame, columns: int=3, spacing: int=20):
    space = "\xa0" * spacing
    k, m = divmod(df.shape[0], columns)
    output = space
    for j in range(columns):
        df_table = df[j*k+min(j, m):(j+1)*k+min(j+1, m)]
        df_styler = df_table.style.set_table_attributes("style='display:inline'")
        output += df_styler._repr_html_() + space
    return output




1.1.4.1 Time taken¶

Measure the time performance of the each scheme and combine into a dataframe.

In [7]:
ls_out = ['']*3
In [8]:
%%timeit -n3
ls_out[0] = simulate_path(ls_scheme['exact'], delta_t, S, arr_rand)[-1]
337 ms ± 3.23 ms per loop (mean ± std. dev. of 7 runs, 3 loops each)
In [9]:
%%timeit -n3
ls_out[1] = simulate_path(ls_scheme['euler_maruyama'], delta_t, S, arr_rand)[-1]
84.1 ms ± 3.43 ms per loop (mean ± std. dev. of 7 runs, 3 loops each)
In [10]:
%%timeit -n3
ls_out[2] = simulate_path(ls_scheme['milstein'], delta_t, S, arr_rand)[-1]
163 ms ± 12.7 ms per loop (mean ± std. dev. of 7 runs, 3 loops each)

From the observed times, it is evident that the Euler-Maruyama scheme is the fastest and is around 4 times faster than the exact solution. While the Milstein scheme took about 50% longer than Euler-Maruyama, it is still a third of the time for the exact solution.


In [11]:
df_sims = pd.DataFrame( ls_out, index=['exact', 'euler_maruyama','milstein']).T
In [12]:
df_sims.describe()
Out[12]:
exact euler_maruyama milstein
count 100000.000000 100000.000000 100000.000000
mean 104.983647 104.984039 104.969628
std 21.154043 21.149681 21.348863
min 41.090237 41.023505 42.191759
25% 89.962204 89.954568 89.803186
50% 102.846399 102.840868 102.722896
75% 117.775404 117.771684 117.738598
max 273.659534 272.468921 290.989523

The different schemes appear to have similar averages but the Milstein scheme has a slightly higher standard deviation and a higher max value.

Notice the average return is around 5% which is what is expected under risk neutral conditions since the replicating portfolio must perform the same as the riskfree asset for the no arbitrage condition to hold. Also, the median is below the mean, therefore expecting right-tailed/right-skewed distributions.


1.1.4.2 Distribution of simulation¶

Examine the histogram of the output.

In [13]:
ls_hist = []
for k in range(df_sims.shape[1]):
    ls_hist.append( np.histogram( df_sims.iloc[:,k], bins=100))
In [14]:
fig1 = go.Figure()
fig1.add_trace(go.Scatter(x=ls_hist[0][1], y=ls_hist[0][0],
                    mode='lines', name='Exact', line=dict(shape='hv'),
                    hovertemplate='Exact - bin: $%{x:.3f} \t count: %{y:,}<extra></extra>'))
fig1.add_trace(go.Scatter(x=ls_hist[1][1], y=ls_hist[1][0],
                    mode='lines', name='E-M', line=dict(shape='hv'),
                    hovertemplate='E-M - bin: $%{x:.3f} \t count: %{y:,}<extra></extra>'))
fig1.add_trace(go.Scatter(x=ls_hist[2][1], y=ls_hist[2][0],
                    mode='lines', name='Milstein', line=dict(shape='hv'),
                    hovertemplate='Milstein - bin: $%{x:.3f} \t count: %{y:,}<extra></extra>'))
fig1.update_yaxes( showspikes=True)    #, range=[104, 105.3] )
fig1.update_xaxes( showspikes=True, tickformat=',')
fig1.update_traces( line_width=3, opacity=0.55)
fig1.update_layout( title='Figure 1 - Histogram of estimated Stock Price', hovermode="x",
                    xaxis_title="Iterations",
                    yaxis_title="Stock price ($)",
                    legend_title="Scheme",
                    width=850, height=550)
fig1.show()

In general, all schemes are right tailed. The output of the Milstein scheme seems to have a longer tail and higher peak compared to the rest. The Euler-Maruyama scheme appears to be tracking the exact solution quite well.

*Note: Easier to observe in the html version as the chart has zoom functionality.

1.1.4.3 Expected price at time $T$¶

Iterate over columns to calculate mean across the number of iterations using expanding mean.

In [15]:
for column in df_sims:
    df_sims[ column + '_mean'] = df_sims[column].expanding().mean()

Look at convergence by taking the expanding mean of dataframe.

In [16]:
fig2 = go.Figure()
fig2.add_trace(go.Scatter(x=df_sims.index, y=df_sims.exact_mean,
                    mode='lines', name='Exact',
                    hovertemplate='Exact: $%{y:.3f}<extra></extra>'))
fig2.add_trace(go.Scatter(x=df_sims.index, y=df_sims.euler_maruyama_mean,
                    mode='lines', name='E-M',
                    hovertemplate='E-M: $%{y:.3f}<extra></extra>'))
fig2.add_trace(go.Scatter(x=df_sims.index, y=df_sims.milstein_mean,
                    mode='lines', name='Milstein',
                    hovertemplate='Milstein: $%{y:.3f}<extra></extra>'))
fig2.update_yaxes( showspikes=True, range=[104, 105.15], tickformat='.1f' )
fig2.update_xaxes( showspikes=True, tickformat=',')
fig2.update_traces( line_width=1, opacity=0.75)
fig2.update_layout( title='Figure 2 - Price vs Iterations', hovermode="x",
                    xaxis_title="Iterations",
                    yaxis_title="Stock price ($)",
                    legend_title="Scheme",
                    width=850, height=550)
fig2.show()

In this case, the difference is very small. However, Euler-Maruyama performs better than the Milstein scheme.

*Note: Chart has been trimmed for better viewing.

1.1.4.4 Error percentage¶

Compare the error of Euler-Maruyama and Milstein to the exact solution.

Define error as the difference between the scheme and the exact, $$ \varepsilon = \hat{S}_{scheme} - S_{exact} $$

In [86]:
df_sims['em_error'] = df_sims.euler_maruyama_mean - df_sims.exact_mean
df_sims['milstein_error'] = df_sims.milstein_mean - df_sims.exact_mean
df_sims['em_error_percentage'] = 100 * (df_sims.em_error / df_sims.exact_mean)
df_sims['milstein_error_percentage'] = 100 * (df_sims.milstein_error / df_sims.exact_mean)
In [87]:
fig3 = go.Figure()
fig3.add_trace(go.Scatter(x=df_sims.index, y=df_sims.em_error_percentage,
                          mode='lines', name='E-M',
                          hovertemplate='<b>E-M error (%):</b> %{y:.4%}<extra></extra>'))
fig3.add_trace(go.Scatter(x=df_sims.index, y=df_sims.milstein_error_percentage,
                          mode='lines', name='Milstein',
                          hovertemplate='<b>Milstein error (%):</b> %{y:.4%}<extra></extra>'))
fig3.update_yaxes( showspikes=True, range=[-0.12, 0.05], tickformat='.1%', dtick=0.02 )
fig3.update_xaxes( showspikes=True)
fig3.update_traces( line_width=2.5, opacity=0.85)
fig3.update_layout( title='Figure 3 - Error comparison', hovermode="x",
                    xaxis_title="Iterations",
                    yaxis_title="Error (%)",
                    legend_title="Scheme",
                    width=850, height=550)
fig3.show()

The error is significantly smaller for Euler-Maruyama while Milstein scheme is still quite good below 2%. However, it is interesting to note that the Milstein scheme does not truly converge to the exact solution. The E-M scheme is within 0.1% after around 7,000 iterations (out of 100,000 iterations).

Based on this, it appears that for this case E-M scheme offers both the speed and accuracy over the Milstein scheme.

*Note: Chart has been trimmed for better viewing.

1.1.4.5 Varying $\delta t$¶

Since the accuracy of the approximation schemes are a function of $\delta t $, consider the effect of varying $\delta t $. This is done by increasing the timesteps (partition size) from daily to $n$ multiple days which gives round (integer) numbers for a year. Increasing step size results in lower number of partitions as the time horizon is maintained to be 1 year.

$$ \text{timestep, } t = \frac{252}{n} \\ $$

For $n=1, \; \delta t = \frac{1}{252}$
For $n=2, \; \delta t = \frac{1}{252 \; / \; 2} = \frac{2}{252}$

and so forth for up to $n=42$ days (two months).

The exact list is given in the table below.

In [88]:
days = 42
n = [ k  for k in range(1,days+1) if 252 % k == 0]
# init dataframe to store results
df_deltaT = pd.DataFrame( np.zeros((len(n),5)),
                          columns=['n', 't','dt','em_error_perc','milstein_error_perc'])
df_deltaT['n'] = n
df_deltaT['t'] = (252/df_deltaT.n).astype('int')
df_deltaT['dt'] = T / df_deltaT.t
df_deltaT.iloc[:,:3]
Out[88]:
n t dt
0 1 252 0.003968
1 2 126 0.007937
2 3 84 0.011905
3 4 63 0.015873
4 6 42 0.023810
5 7 36 0.027778
6 9 28 0.035714
7 12 21 0.047619
8 14 18 0.055556
9 18 14 0.071429
10 21 12 0.083333
11 28 9 0.111111
12 36 7 0.142857
13 42 6 0.166667
In [89]:
# init empty stock price array - init to largest required size, then slice as necessary for larger timesteps
S_dt =  np.zeros(( max(df_deltaT.t), runs))
S_dt[0] = S0

# initialize random values
np.random.seed(20230419)
arr_rand_dt = np.random.standard_normal(max(df_deltaT.t)*runs).reshape( max(df_deltaT.t),runs)
In [90]:
for k in df_deltaT.index:
    df_dt = pd.DataFrame({
            'exact': simulate_path( ls_scheme['exact'], df_deltaT.dt[k], S_dt[:df_deltaT.t[k]],
                                    arr_rand_dt[:df_deltaT.t[k]])[-1],
            'euler_maruyama': simulate_path( ls_scheme['euler_maruyama'], df_deltaT.dt[k],
                                             S_dt[:df_deltaT.t[k]], arr_rand_dt[:df_deltaT.t[k]])[-1],
            'milstein': simulate_path( ls_scheme['milstein'], df_deltaT.dt[k], S_dt[:df_deltaT.t[k]],
                                       arr_rand_dt[:df_deltaT.t[k]])[-1],
    })
    df_deltaT.loc[k,'em_error_perc'] = (100 * (np.mean(df_dt.euler_maruyama) - np.mean(df_dt.exact)) 
                                        / np.mean(df_dt.exact)).item()
    df_deltaT.loc[k,'milstein_error_perc'] = (100 * (np.mean(df_dt.milstein) - np.mean(df_dt.exact)) 
                                              / np.mean(df_dt.exact)).item()
    
In [91]:
display_html( style_table( df_deltaT, 2), raw=True)
                    
  n t dt em_error_perc milstein_error_perc
0 1 252 0.003968 0.000373 -0.013354
1 2 126 0.007937 0.000455 -0.015436
2 3 84 0.011905 -0.000245 -0.011239
3 4 63 0.015873 -0.000210 -0.013360
4 6 42 0.023810 -0.002310 -0.005157
5 7 36 0.027778 -0.002745 -0.005314
6 9 28 0.035714 -0.002499 -0.011095
                    
  n t dt em_error_perc milstein_error_perc
7 12 21 0.047619 -0.004517 -0.008735
8 14 18 0.055556 -0.006080 -0.007073
9 18 14 0.071429 -0.008485 -0.006405
10 21 12 0.083333 -0.010361 -0.006508
11 28 9 0.111111 -0.012498 -0.011254
12 36 7 0.142857 -0.017226 -0.010987
13 42 6 0.166667 -0.017955 -0.015620
                    

Table above showing the varying $\delta t$ and error of schemes.



In [92]:
fig4 = go.Figure()
fig4.add_trace(go.Scatter(x=df_deltaT.dt, y=df_deltaT.em_error_perc,
                          mode='lines+markers', name='E-M', 
                          hovertemplate='<b>dt:</b> %{x:.4f}<br>' + f'<b>E-M error (%):</b> ' +
                                          '%{y:.4%}<extra></extra>'
                         ))
fig4.add_trace(go.Scatter(x=df_deltaT.dt, y=df_deltaT.milstein_error_perc, 
                          mode='lines+markers', name='Milstein', 
                          hovertemplate='<b>dt:</b> %{x:.4f}<br>' + f'<b>Milstein error (%):</b> ' +
                                          '%{y:.4%}<extra></extra>'
                         ))
fig4.update_yaxes( showspikes=True, tickformat='.1%', dtick=0.002 ) #range=[-0.55, 0.015], 
fig4.update_xaxes( showspikes=True, tickformat='0.2f')
fig4.update_layout( title='Figure 4 - Error (%) vs Iterations' + 
                          '<br><sup>for varying dt at 100,000 iterations</sup>',
                    hovermode="x", xaxis_title="dt", yaxis_title="Error (%)",
                    legend_title="Scheme", width=850, height=550)
fig4.show()

It is expected that as $\delta t$ increases, the error would increase. The error for Euler-Maruyama scheme appears to grow somewhat linearly. While the error of the Milstein scheme seems to be improving initially but as the step size increases, the error starts to increase as well. This error turn around is most likely driven by the correction term which has a $\sigma^2$. More importantly, these errors are at the second/third decimal places onwards which is already very small. In most cases, any decision will most likely be driven by the computing resources to run these calculations.

*Note: Negative sign indicates underestimating.


1.1.5 Conclusion for GBM schemes¶

In conclusion, from all the comparisons above, it appears that the preferred scheme will be between the exact solution and the Euler-Maruyama scheme depending on speed, low complexity and low error. For small computations, the exact solution will still be performant, however as the simulation size grows there will be a need to improve performance.



1.2 Asian Option Value¶

For an Asian option, the expected value of the discounted payoff under the risk neutral density $\mathbb{Q}$ is given by,

$$ V (S, t) = e^{-r(T-t)} \; \mathbb{E^\mathbb{Q}}[ Payoff] $$

For call and put options, consider the following:

  1. Fixed vs floating strike
  2. Discrete vs continuous sampling
  3. Arithmetic vs geometric averaging

Using the previously provided data, $$ \begin{align} \text{Today’s stock price, } S_0 &= 100 \\[5pt] \text{Time, } T &= 1 \text{ year} \\[5pt] \text{volatility, } \sigma &= 20\% \\[5pt] \text{constant risk-free interest rate, } r &= 5\% \\[5pt] \end{align} $$

and

$$ \begin{align} \text{Strike, } E &= 100 \\[5pt] \text{Time to expiry, } (T-t) &= 1 \text{ year} \\[5pt] \end{align} $$

It is important to understand that the stock price is a Martingale under the $\mathbb{Q}$ measure and the mean or drift rate is zero. This is similar to what was considered in section 1.1.

For a floating strike, the strike price is not pre-determined but instead is calculated as a function of the stock price which means it will move as the stock price moves. Clearly, this makes the option heavily path dependent as the path the stock price takes will determine the payoff.

The calculation of the strike price may take almost any form such as weighted, exponential weighted, volume weighted and so forth. In this case this analysis will compare between geometric and arithmetic means.


1.2.1 Payoff Function¶

In general, the payoff function for an option has the form,

$$ \begin{align} P_{call} &= max(S - E, \; 0) \\[4pt] \end{align} $$

and a put option,

$$ \begin{align} P_{put} &= max(E - S, \; 0) \end{align} $$ $$ \text{where } \quad S=\text{Stock price at expiry} \quad \text{and} \quad E = \text{strike price} \\ $$

These are the exact payoff function for a fixed strike call and put option (European options). None of these are path dependent as the payoff only depends on the price of the underlying asset at expiry.


For options with a floating strike, whether arithmetic or geometric, continuous or discrete payoff, it will involve modifying the strike price, $E$.

First, denote the floating strike as the function $I(S)$, hence the payoff function of the form,

$$ P_{call, floating} = max(S - I, 0) \\[5pt] $$

and a put option,

$$ P_{put, floating} = max(I - S, 0) \\[5pt] $$

and look at the different forms of $I(S)$.


Floating arithmetic payoff $$ I_{arithmetic}(S_t) \;=\; \frac{1}{T} \int_{0}^{T} S(\tau) \; d\tau \;=\; \frac{1}{N} \sum_{i=0}^{N} S_i \\[5pt] $$


Floating geometric payoff $$ I_{geometric}(S) = \left( \prod_{i = 0}^{N} S_i \right)^\frac{1}{N} \;=\; \sqrt[N]{S_0 * S_1*...*S_N} \;=\; exp \left\{ \frac{1}{N} \sum_{i=0}^{N} \; ln \; S_i \right\} \\[5pt] $$

For continuously sampled, $N$ is the number of days to expiry whereas for the discrete case, $N$ is the number of dates sampled.

While in theory it is possible to consider the stock price in continuous time, in the real world, this is not possible. So the stock price is discretized by sampling and this may take the form of continuous sampling at some constant rate (such as daily) or discrete sampling where the stock price is sampled on specific dates or larger than daily constant interval (such as weekly or monthly).

Of course, each of these affects the strike price in different ways and more importantly, maps how the path of the underlying asset will affect the strike price and in turn, the payoff at expiry.




In the following analysis,

  1. the discrete case will consider the sampling at weekly ($n=5$), fortnightly ($n=10$) and monthly ($n=21$).
  2. while the continuous will sample at the daily level.
  3. Use of exact solution for stock prices.

1.2.2 Simulation for call and put option¶

Use the previous defined functions and variables, begin by defining the array of underlying prices, the discount term and a function to calculate value of the option.

In [93]:
# Monte Carlo Price
arr_underlying = simulate_path(ls_scheme['exact'], delta_t, S, arr_rand)

strike = 100.0
tau = 1

discount_term = np.exp( -rate*tau )
dict_discrete_sampling = {'weekly': 5, 'fortnightly': 10, 'monthly': 21}

# store results for plotting
df_call = pd.DataFrame({'type':[], 'payoff': [], 'value': []})
df_put = pd.DataFrame({'type':[], 'payoff': [], 'value': []})
dict_payoff = { 'call': {}, 'put': {}, 'hist_call': {}, 'hist_put': {} }

Define functions to perform strike calculations, expected payoff and expected value of discounted payoff.

In [94]:
def fixed_strike(strike, size) -> np.array:
    return np.repeat(strike, size)

def arith_cont(arr_price: np.array) -> np.array:
    return np.cumsum(arr_price, axis=0) / np.arange(1, arr_price.shape[0]+1)[:,np.newaxis]

def geom_cont(arr_price: np.array) -> np.array:
    return np.exp(( np.cumsum( np.log(arr_price), axis=0) / 
                    np.arange(1, arr_price.shape[0]+1)[:,np.newaxis])[-1])

def arith_discrete(arr_price: np.array, n: int) -> np.array:
    return ( np.cumsum(arr_price[::n], axis=0) / 
             np.arange(1, arr_price[::n].shape[0]+1)[:,np.newaxis])

def geom_discrete(arr_price: np.array, n: int) -> np.array:
    return np.exp(( np.cumsum( np.log(arr_price[::n]), axis=0) /
                    np.arange(1, arr_price[::n].shape[0]+1)[:,np.newaxis])[-1])

# expected payoff and expected value of discounted payoff
def value_option( arr_price: np.array, arr_strike: float, discounting: float, option_type: str="" ) -> np.array:
    if option_type.lower() == 'call':
        arr_payoff = arr_price[-1] - arr_strike
    elif option_type.lower() == 'put':
        arr_payoff = arr_strike - arr_price[-1]
    arr_payoff_filt = np.ma.filled(np.ma.masked_less_equal(arr_payoff, 0), fill_value=0)
    option_payoff = np.mean( arr_payoff_filt)
    option_value = discounting * option_payoff
    return option_payoff, option_value

# calculate all strikes
def calc_strikes( strike, arr_price: np.array, discrete_sample: dict) -> dict:
    dict_strikes = { 'fixed': fixed_strike(strike, arr_underlying.shape[1]),
                    'arithmetic continuous': arith_cont(arr_underlying)[-1],
                    'geom continuous': geom_cont(arr_underlying) }
    for key,value in dict_discrete_sampling.items():
        dict_strikes[f'arithmetic discrete {key}'] = arith_discrete(arr_underlying, value)[-1]
        dict_strikes[f'geom discrete {key}'] = geom_discrete(arr_underlying, value)
    return dict_strikes
In [95]:
dict_strike = calc_strikes( strike, arr_underlying, dict_discrete_sampling)

# calculate expected payoff and expected value of discounted payoff
for i,key in enumerate(dict_strike.keys()):
    df_call.loc[i,:] = (key,) + value_option( arr_underlying, dict_strike[key], discount_term, 'call')
    df_put.loc[i,:] = (key,) + value_option( arr_underlying, dict_strike[key], discount_term, 'put')
In [96]:
dfcall_styler = df_call.sort_values('value', ascending=False
                      ).style.set_table_attributes("style='display:inline'"
                      ).set_caption('Call Options for various strike')
dfput_styler = df_put.sort_values('value', ascending=False
                    ).style.set_table_attributes("style='display:inline'"
                    ).set_caption('Put Options for various strike')
display_html( "\xa0"*20 + dfcall_styler._repr_html_() + "\xa0" * 30 + dfput_styler._repr_html_(), raw=True)
                    
Call Options for various strike
  type payoff value
0 fixed 10.876316 10.345871
8 geom discrete monthly 6.707436 6.380310
7 arithmetic discrete monthly 6.492180 6.175552
4 geom discrete weekly 6.327508 6.018911
2 geom continuous 6.319354 6.011156
6 geom discrete fortnightly 6.308866 6.001179
3 arithmetic discrete weekly 6.106878 5.809042
1 arithmetic continuous 6.101603 5.804024
5 arithmetic discrete fortnightly 6.083495 5.786799
                              
Put Options for various strike
  type payoff value
0 fixed 5.892668 5.605279
7 arithmetic discrete monthly 3.762504 3.579004
8 geom discrete monthly 3.630305 3.453253
1 arithmetic continuous 3.571504 3.397319
3 arithmetic discrete weekly 3.567416 3.393431
5 arithmetic discrete fortnightly 3.544841 3.371957
2 geom continuous 3.438694 3.270987
4 geom discrete weekly 3.433377 3.265929
6 geom discrete fortnightly 3.408852 3.242600

Table above displays the expected payoff and expected value of the discounted payoff for both call and put options for the different types of strikes ordered by descending value.


1.2.3 Various Strike Call Option¶

In [97]:
fig5 = go.Figure()
for columns in df_call.T:
    fig5.add_trace(go.Scatter(x=[0,1], y=df_call.T.iloc[:0:-1,columns], mode='lines+markers', 
                              name=f'{df_call.T.iloc[0,columns]}', meta=f'{df_call.T.iloc[0,columns]}',
                              hovertemplate='<b>%{meta}</b> - $%{y:.3f}<extra></extra>',
                              line=dict(width=1.3, dash='solid'),
                              marker=dict(size=8, symbol='x')
                             ))
fig5.update_yaxes( showspikes=True)
fig5.update_xaxes( showspikes=True)
fig5.update_layout( title='Figure 5 - Monte Carlo <b>Call</b> Option Estimation' +
                          '<br><sup>100,000 iterations - various strike types</sup>',
                    hovermode="x", xaxis_title="Time", yaxis_title="Payoff ($)",
                    legend_title="Strike type", width=850, height=475)
fig5.show()
The chart above displays the expected payoff and the expected value of the discounted payoff at $t=0$ and $t=T=1$. The relationship is simply the discounting to the starting time.

Among the Pythagorean means, the arithmetic mean is the largest followed by the geometric mean. Therefore, it is expected that the options with arithmetic averaging will have a lower payoff (higher strike) compared to geometric averaging. As expected, the fixed strike option will have the highest value. The limited data for discrete arithmetic mean results in a lower strike and in turn higher payoff as compared to the continuous arithmetic mean. When comparing discrete sampling regimes, monthly sampling has higher payoff versus fornightly and weekly. Essentially, the less data used for sampling results in lower strike and higher payoff. In effect allowing volatility to drive prices higher than the sampled average. This line of thinking could easily be extended to include the fixed strike as well where it has no sampling.

For the fixed strike, there is no path dependence, therefore the strike is the level for which the underlying either exceeds or does not reach at expiry. Where as for the rest, the path at the sample points dictates the strike.

Comparisons which results in higher call option value in order of significance:

  1. fixed > floating
  2. low sampling > higher sampling
  3. geometric > arithmetic

1.2.4 Various Strike Put Option¶

In [98]:
fig7 = go.Figure()
for columns in df_put.T:
    fig7.add_trace(go.Scatter(x=[0,1], y=df_put.T.iloc[:0:-1,columns], mode='lines+markers', 
                              name=f'{df_put.T.iloc[0,columns]}', meta=f'{df_put.T.iloc[0,columns]}',
                              hovertemplate='<b>%{meta}</b> - $%{y:.3f}<extra></extra>',
                              line=dict( width=1.3, dash='solid'),
                              marker=dict(size=8, symbol='x')
                             ))
fig7.update_yaxes( showspikes=True)
fig7.update_xaxes( showspikes=True)
fig7.update_layout( title='Figure 7 - Monte Carlo <b>Put</b> Option Estimation' +
                          '<br><sup>100,000 iterations - various strike types</sup>',
                    hovermode="x", xaxis_title="Time", yaxis_title="Payoff ($)",
                    legend_title="Strike type", width=870, height=475)
fig7.show()

The chart above shows the expected payoff and the expected value of the discounted payoff at $t=0$ and $t=T=1$. Similar to the chart for call options, the relationship is defined by the discounting factor.

The first thing that is noticeable is that the values for the floating strike are approximately half of the fixed strike. The fixed strike has the highest value.

In the case of the put option, it is desirable to have strike as large as possible since the payoff function is, $P = I - S$. Therefore for puts, arithmetic mean results in higher value than geometric averaging.

For sampling time, it is similar to call options where less sampling results in higher value.

Comparisons which results in higher put option value in order of significance:

  1. fixed > floating
  2. low sampling > higher sampling
  3. arithmetic > geometric

1.2.5 Varying values¶

Previously, the effect of varying $\delta t$ has been investigated. This section will vary the data to discover the effect on option price. Specifically the volatility and time to expiry. The expected effects are:

  1. volatility - more movement in the underlying price may result in higher likelihood of being in the money
  2. time to expiry - longer time to expiry possibly results in more opportunity for the underlying to move

1.2.5.1 Varying Volatility¶

In [99]:
ls_sigma = np.linspace(0.02,0.4,20) #np.linspace(0.1,0.6,11)
colnames_vol = ['vols'] + list(dict_strike.keys())
dict_vols = { 'call': pd.DataFrame( columns=colnames_vol),
              'put': pd.DataFrame( columns=colnames_vol) }
dict_vols['call'].loc[:,'vols'] = dict_vols['put'].loc[:,'vols'] = ls_sigma
In [100]:
ls_sigma
Out[100]:
array([0.02, 0.04, 0.06, 0.08, 0.1 , 0.12, 0.14, 0.16, 0.18, 0.2 , 0.22,
       0.24, 0.26, 0.28, 0.3 , 0.32, 0.34, 0.36, 0.38, 0.4 ])
In [101]:
for j,vol in enumerate(ls_sigma):
    arr_underlying1 = simulate_path( partial(exact, r=rate, vol=vol), delta_t, S, arr_rand)
    dict_strike = calc_strikes( strike, arr_underlying1, dict_discrete_sampling)
    for i,key in enumerate(dict_strike.keys()):
        dict_vols['call'].loc[j, key] = value_option( arr_underlying1, dict_strike[key], 
                                                      discount_term, 'call')[1]
        dict_vols['put'].loc[j, key] = value_option( arr_underlying1, dict_strike[key], 
                                                     discount_term, 'put')[1]
In [102]:
fig8 = go.Figure()
for column in dict_vols['call'].drop('vols', axis=1):
    fig8.add_trace(go.Scatter(x=dict_vols['call'].vols, y=dict_vols['call'].loc[:,column], mode='lines+markers', 
                              name=f'{column}', meta=f'{column}',
                              hovertemplate='%{meta}: $%{y:.3f}<extra></extra>',
                              line=dict( width=1.3, dash='solid'),
                              marker=dict(size=8, symbol='x')
                             ))
fig8.update_yaxes( showspikes=True)
fig8.update_xaxes( showspikes=True, dtick=0.04, tickformat='.2f')
fig8.update_layout( title='Figure 8 - Call Option Value for Varying volatility' +
                          '<br><sup>100,000 iterations - various strike types</sup>',
                    hovermode="x", xaxis_title="volatility", yaxis_title="Value ($)",
                    legend_title="Strike type", width=870, height=600)
fig8.show()

In general, it is expected that as volatility increases, the option value increases as well. Since low volatility indicates that the underlying is less likely to move whereas a higher volatility means that the underlying will be move higher or lower at each step. Interestingly, for the floating strike options there seems to be some second order relationship at volatilities below 0.2 and beyond 0.2, the relationship tends to be linear. This behaved is observed in the fixed strike as well although much less pronounced. Potentially this is due to the riskfree term $r - \frac{1}{2} \sigma^2$ in the underlying and as $\sigma$ becomes larger, the diffusion term takes over. The floating strike in a way amplifies this as the averaging is able to keep up with volatility is low.



In [103]:
fig9 = go.Figure()
for column in dict_vols['put'].drop('vols', axis=1):
    fig9.add_trace(go.Scatter(x=dict_vols['put'].vols, y=dict_vols['put'].loc[:,column], mode='lines+markers', 
                              name=f'{column}',
                              meta=f'{column}',
                              hovertemplate='%{meta}: $%{y:.3f}<extra></extra>',
                              line=dict( width=1.3, dash='solid'),
                              marker=dict(size=8, symbol='x')
                             ))
fig9.update_yaxes( showspikes=True)
fig9.update_xaxes( showspikes=True, dtick=0.04, tickformat='.2f')
fig9.update_layout( title='Figure 9 - Put Option Value for Varying volatility' +
                          '<br><sup>100,000 iterations - various strike types</sup>',
                    hovermode="x", xaxis_title="volatility", yaxis_title="Value ($)",
                    legend_title="Strike type", width=870, height=600)
fig9.show()

Similar to the call option, increasing volatility increases the put option value. However, it is interesting to note that the option value is almost zero as the underlying price will not have move significantly. Yet, for the Asian option floating strike, the value appears to be higher for low volatility and has the same second order hockey stick profile.



1.2.5.2 Varying Time to Expiry¶

Investigate the effect of varying time to expiry while keeping all other data constant, specifically $\delta t$ shall remain at 1 day. The time to expiry being considered as in increments of 2 months from 1 months up till 18 months.

In [104]:
months = [ i for i in range(1, 19, 1)] 

colnames_mte = ['months', 'horizon', 'step', 'discount'] + list(dict_strike.keys())
dict_mte = { 'call': pd.DataFrame( columns=colnames_mte),
              'put': pd.DataFrame( columns=colnames_mte)}
for key in dict_vols.keys():
    dict_mte[key].loc[:,'months'] = months
    dict_mte[key].loc[:,'horizon'] = dict_mte[key].months / 12
    dict_mte[key].loc[:,'step'] = (dict_mte[key].horizon * 252).astype('int')
    dict_mte[key].loc[:,'discount'] = np.exp( (-rate*dict_mte[key].horizon).to_list())
In [105]:
dict_mte['call'].iloc[:,:4]
Out[105]:
months horizon step discount
0 1 0.083333 21 0.995842
1 2 0.166667 42 0.991701
2 3 0.25 63 0.987578
3 4 0.333333 84 0.983471
4 5 0.416667 105 0.979382
5 6 0.5 126 0.97531
6 7 0.583333 147 0.971255
7 8 0.666667 168 0.967216
8 9 0.75 189 0.963194
9 10 0.833333 210 0.959189
10 11 0.916667 231 0.955201
11 12 1.0 252 0.951229
12 13 1.083333 273 0.947274
13 14 1.166667 294 0.943335
14 15 1.25 315 0.939413
15 16 1.333333 336 0.935507
16 17 1.416667 357 0.931617
17 18 1.5 378 0.927743
In [106]:
# initialize arrays to largest required size, then slice as necessary for larger timesteps

# initialize empty stock price array 
S_mte =  np.zeros(( max(dict_mte['call'].step), runs))
S_mte[0] = S0

# initialize random values
np.random.seed(20230419)
arr_rand_mte = np.random.standard_normal(max(dict_mte['call'].step)*runs
                                        ).reshape( max(dict_mte['call'].step),runs)
In [107]:
for j in dict_mte['call'].index:
    arr_underlying2 = simulate_path( partial(exact, r=rate, vol=sigma), delta_t,
                                     S_mte[:dict_mte['call'].step[j]],
                                     arr_rand_mte[:dict_mte['call'].step[j]])
    dict_strike = calc_strikes( strike, arr_underlying2, dict_discrete_sampling)
    for i,key in enumerate(dict_strike.keys()):
        dict_mte['call'].loc[j, key] = value_option( arr_underlying2, dict_strike[key], 
                                                      dict_mte['call'].discount[j], 'call')[1]
        dict_mte['put'].loc[j, key] = value_option( arr_underlying2, dict_strike[key], 
                                                     dict_mte['call'].discount[j], 'put')[1]
In [108]:
fig10 = go.Figure()
for column in dict_mte['call'].iloc[:,4:]:
    fig10.add_trace(go.Scatter(x=dict_mte['call'].months, y=dict_mte['call'].loc[:,column],
                               mode='lines+markers', name=f'{column}', meta=f'{column}',
                               hovertemplate='%{meta}: $%{y:.3f}<extra></extra>',
                               line=dict( width=1.3, dash='solid'),
                               marker=dict(size=8, symbol='x')
                             ))
fig10.update_yaxes( showspikes=True)
fig10.update_xaxes( showspikes=True, dtick = 3)
fig10.update_layout( title='Figure 10 - Call Option Value for Varying Time to Expiry' +
                          '<br><sup>100,000 iterations - various strike types</sup>',
                    hovermode="x", xaxis_title="months", yaxis_title="Value ($)",
                    legend_title="Strike type", width=870, height=600)
fig10.show()

As expected, time to expiry is proportional to the value of the option. Given that the underlying has sufficiently time to diffuse or move through its price range. Option value is also influenced by the discounting factor relationship $e^{-rT}$ as the time to expiry increases. In the chart above, the hockey stick profile is also visible


In [109]:
fig11 = go.Figure()
for column in dict_mte['put'].iloc[:,4:]:
    fig11.add_trace(go.Scatter(x=dict_mte['put'].months, y=dict_mte['put'].loc[:,column],
                               mode='lines+markers', name=f'{column}', meta=f'{column}',
                               hovertemplate='%{meta}: $%{y:.3f}<extra></extra>',
                               line=dict( width=1.3, dash='solid'),
                               marker=dict(size=8, symbol='x')
                             ))
fig11.update_yaxes( showspikes=True)
fig11.update_xaxes( showspikes=True, dtick = 3)
fig11.update_layout( title='Figure 11 - Put Option Value for Varying  Time to Expiry' +
                          '<br><sup>100,000 iterations - various strike types</sup>',
                    hovermode="x", xaxis_title="months", yaxis_title="Value ($)",
                    legend_title="Strike type", width=870, height=600)
fig11.show()

The chart above for put options indicates that for short expiry options (2 months and below) for this set of data, the floating strike scheme has high option value with high sampling arithmetic averaging being highest. Again, this is potentially due to the low data samples for averaging.



1.2.6 Conclusion¶

The use of closed form solution, Euler-Maruyama or Milstein scheme will depend on the use case. However, for running Monte Carlo simulation of stock prices, the Milstein scheme does not seem to offer any advantage as the Euler-Maruyama appears to offer the trifecta of speed, better accuracy and faster convergence.

The comparisons above highlight the different methods for strike calculation and it is clear in all cases that a fixed strike (European option), with all else being equal, should be preferred as there are no path dependence. However, as demonstrated, path dependence can be taken into account in the simulation.

In order of highest value for options:

  1. fixed > floating
  2. high volatility > low volatility
  3. longer duration
  4. low sampling > higher sampling
  5. geometric > arithmetic (for call options and the reverse is true for put options)

Further investigation could be done for floating price options as well as other method for floating strike calculations such as volume weighted average or exponential weighted average.





Question 2¶

The model problem is given as

$$ \frac{d^2y}{dx^2} + P\;(x) \frac{dy}{dx} + Q\,(x)\;y = f\,(x) \\[15pt] \text{subject to the boundary conditions,} \quad y\;(a) = \alpha \quad \text{and} \quad y\;(b) = \beta \\ $$

and $x$ is given as a regular partition of the interval $[a,b]$ where,

$$ x_0 < x_1 < x_2 <....< x_{n-1} < x_n \\[5pt] \text{such that} \\ x_0 = a \quad \text{and} \quad x_n = b \\ $$

with the partition size,

$$ \delta x = \frac{b-a}{n} $$

Note that the model has a single input variable, $x$.



2.1 Numerical Approximation¶

Let, $$ y_i \;=\; y\;(x_i) \\ P_i \;=\; P\;(x_i) \\ Q_i \;=\; Q\;(x_i) \\ f_i \;=\; f\;(x_i) \\ $$

Using central difference approximation, start by applying Taylor Series Expansion on $y\,(x)$, for $ x + \delta x $ around $ x $

$$ y\;(x + \delta x) \;=\; y\;(x) \;+\; \delta x \frac{dy}{dx} \;+\; \frac{\delta x^2}{2} \frac{d^2y}{dx^2} \;+\; \frac{\delta x^3}{3!} \frac{ d^3 y}{ dx^3} \;+\; O\;(\delta x^4) \tag{1} \\ $$

and similarly for $ x - \delta x $ around $ x $

$$ y\;(x - \delta x) \;=\; y\;(x) \; - \; \delta x \frac{dy}{dx} \;+\; \frac{\delta x^2}{2} \frac{d^2y}{dx^2} \;-\; \frac{\delta x^3}{3!} \frac{ d^3 y}{ dx^3} \;+\; O\;(\delta x^4) \tag{2} \\ $$

Consider the difference of $\;(1) \;-\; (2)$,

$$ \begin{align} y\;(x+ \delta x) \;-\; y\;(x-\delta x) \;&=\; 2\; \delta x\; \frac{dy}{dx} + O\;(\delta x^3) \\[8pt] \frac{dy}{dx} \;&=\; \frac{y\;(x+ \delta x) \;-\; y\;(x-\delta x)}{ 2\; \delta x\; } \\[8pt] \;&=\; \frac{ y_{i+1} \;-\; y_{i-1} }{ 2 \; \delta x } \qquad for \quad i = 1,2,...,n-1\\[8pt] \end{align} $$

Consider the sum of $\;(1) \;+\; (2)$,

$$ \begin{align} y\;(x+ \delta x) \;+\; y\;(x-\delta x) \;&=\; 2y\;(x) \;+\; \delta x^2 \; \frac{d^2 y}{dx^2} \;+\; O\;(\delta x^4) \\[8pt] \frac{ y\;(x+ \delta x) \;-\; 2y\;(x) \;+\; y\;(x-\delta x) } { \delta x^2 } \;&=\; \frac{ d^2 y }{ dx^2 } \;+\; O\;(\delta x^4) \\[8pt] \frac{d^2 y}{dx^2} \;&=\; \frac{ y\;(x+ \delta x) \;-\; 2y\;(x) \;+\; y\;(x-\delta x) }{ \delta x^2 } \\[8pt] \;&=\; \frac{ y_{i+1} \;-\; 2y_i \;+\; y_{i-1} }{ \delta x^2 } \qquad for \quad i = 1,2,...,n-1\\[8pt] \end{align} $$



Dr Ahmad, R. (2023). Module 3 Lecture 4 - "Introduction to Numerical Methods"[PowerPoint slides].
CQF Program January 2023



Next, substitute the numerical approximations of the derivatives found above into the model.

$$ \frac{ y_{i+1} \;-\; 2y_i \;+\; y_{i-1} }{ \delta x^2 } \;+\; P_i \left( \; \frac{ y_{i+1} \;-\; y_{i-1} }{ 2 \; \delta x } \; \right) \;+\; Q_i y_i = f_i \\[10pt] \begin{align} y_{i+1} \;-\; 2y_i \;+\; y_{i-1} \;+\; \frac{ \delta x }{2} P_i y_{i+1} \; \;-\; \frac{ \delta x }{2} P_i y_{i-1} \; \;+\; \delta x^2 Q_i y_i \;&=\; \delta x^2 f_i \\[8pt] y_{i+1} \left( 1 \;+\; \frac{ \delta x }{2} P_i \right) \;+\; y_i \left( \;-2 \;+\; \delta x^2 Q_i y_i \right) \;+\; y_{i-1} \left( 1 \;-\; \frac{ \delta x }{2} P_i \right) \;&=\; \delta x^2 f_i \tag{3} \\[15pt] \end{align} $$$$ for \quad i = 1,2,...,n-1 $$




Write the coefficients as,

$$ A_i = \left( 1 + \frac{ \delta x }{2} P_i \right) \\[12pt] B_i = \left( -2 + \delta x^2 Q_i \right) \\[8pt] C_i = \left( 1 - \frac{ \delta x }{2} P_i \right) \\[8pt] $$



2.2 Matrix Representation¶

Equation 3 could be represented in matrix form of linear systems

From the given boundary conditions of $$ y_0 \;=\; y\;(x_0) \;=\; y\;(a) \;=\; \alpha \\ y_n \;=\; y\;(x_n) \;=\; y\;(b) \;=\; \beta \\ $$

Let,

$$ x \; = \; \begin{bmatrix} y\;(x_0) \\ y\;(x_1) \\ y\;(x_2) \\ . \\ . \\ . \\ y\;(x_{n-1}) \\ y\;(x_n) \end{bmatrix} \; = \; \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ . \\ . \\ . \\ y_{n-1} \\ y_n \end{bmatrix} \qquad\qquad b \text{ } = \text{ } \delta x^2 \text{ } \begin{bmatrix} f_0 \\ f_1 \\ f_2 \\ . \\ . \\ . \\ f_{n-1} \\ f_n \end{bmatrix} \; = \; \quad \delta x^2 \begin{bmatrix} \alpha \; \delta x^{-2} \\ f_1 \\ f_2 \\ . \\ . \\ . \\ f_{n-1} \\ \beta \; \delta x^{-2} \end{bmatrix} \\ $$

The model expressed as $ A \boldsymbol{x} = \boldsymbol{b} $ for any arbitrary value of $n$ with the boundary conditions in the first and last row,

$$ \begin{align} \\ \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & . & . & 0 & 0 & 0 & 0 \\ C_1 & B_1 & A_1 & . & . & . & . & . & . & . & . \\ 0 & C_1 & B_1 & A_1 & . & . & . & . & . & . & . \\ . & 0 & C_1 & B_1 & A_1 & . & . & . & . & . & . \\ . & . & 0 & 0 & 0 & . & . & . & . & . & . \\ . & . & . & . & . & . & . & . & . & . & . \\ . & . & . & . & . & . & . & C_{n-2} & B_{n-2} & A_{n-2} & 0 & \\ . & . & . & . & . & . & . & 0 & C_{n-1} & B_{n-1} & A_{n-1} \\ . & . & . & . & . & . & . & . & 0 & 0 & 1 \\ \end{bmatrix} \quad \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ . \\ . \\ . \\ . \\ y_{n-1} \\ y_n \end{bmatrix} \quad &= \quad \delta x^2 \begin{bmatrix} \alpha \; \delta x^{-2} \\ f_1 \\ f_2 \\ . \\ . \\ . \\ . \\ f_{n-1} \\ \beta \; \delta x^{-2} \end{bmatrix} \\ \text{ } \\ A \qquad\qquad\qquad\qquad\qquad \boldsymbol{x} \qquad &= \qquad \boldsymbol{b} \end{align} $$

This could be solved as a matrix inversion problem for $\boldsymbol{x}$ to yield,

$$ \boldsymbol{x} \;=\; A^{-1} \boldsymbol{b} $$



Dr Ahmad, R. (2023). Module 3 Lecture 7 - "Further Finite Difference Methods"[PowerPoint slides].
CQF Program January 2023



2.3 Numerical calculations¶

Given the differential equation,

$$ \frac{d^2y}{dx^2} + 3 \frac{dy}{dx} + 2y = 4x^2 \\[15pt] \text{subject to the boundary conditions,} \quad y\;(1) = 1 \quad \text{and} \quad y\;(2) = 6 \\ $$

Comparing to the model

$$ \frac{d^2y}{dx^2} + P\;(x) \frac{dy}{dx} + Q\,(x)\;y = f\,(x) \\[15pt] \text{subject to the boundary conditions,} \quad y\;(a) = \alpha \quad \text{and} \quad y\;(b) = \beta \\ $$

Obtain values as,

$$ \begin{align} P(x) &= 3 \\ Q(x) &= 2 \\ f(x) &= 4x^2 \\ \end{align} $$$$ a = 1 \quad\Rightarrow\quad y_0 = \alpha = 1 \\ b = 2 \quad\Rightarrow\quad y_n = \beta = 6 \\ $$

Write the coefficients as,

$$ A_i \;=\; 1 + \frac{ \delta x }{2} P_i \quad =\; 1 + \frac{ \delta x }{2} \;(3) \\[12pt] B_i \;=\; -2 + \delta x^2 Q_i \quad =\; -2 + \delta x^2 \; (2) \\[8pt] C_i \;=\; 1 - \frac{ \delta x }{2} P_i \quad =\; 1 - \frac{ \delta x }{2} \; (3) \\[8pt] $$

The matrix equation becomes,

$$ \boldsymbol{x} \; = \; \begin{bmatrix} y\;(x_0) \\ y\;(x_1) \\ y\;(x_2) \\ . \\ . \\ . \\ y\;(x_{n-1}) \\ y\;(x_n) \end{bmatrix} \; = \; \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ . \\ . \\ . \\ y_{n-1} \\ y_n \end{bmatrix} $$
$$ \boldsymbol{b} \; = \; \delta x^2 \begin{bmatrix} f_0 \\ f_1 \\ f_2 \\ . \\ . \\ . \\ f_{n-1} \\ f_n \end{bmatrix} \; = \; \begin{bmatrix} \alpha \\ \delta x^2 f(x_1) \\ \delta x^2 f(x_2) \\ . \\ . \\ . \\ \delta x^2 f(x_{n-1}) \\ \beta \end{bmatrix} \; = \; \begin{bmatrix} 1 \\ 4x_1^2 \delta x^2 \\ 4x_2^2 \delta x^2 \\ . \\ . \\ . \\ 4x_{n-1}^2 \delta x^2 \\ 6 \end{bmatrix} \\ $$
$$ \begin{align} \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & . & . & 0 & 0 & 0 & 0 \\ C_1 & B_1 & A_1 & . & . & . & . & . & . & . & . \\ 0 & C_1 & B_1 & A_1 & . & . & . & . & . & . & . \\ . & 0 & C_1 & B_1 & A_1 & . & . & . & . & . & . \\ . & . & 0 & 0 & 0 & . & . & . & . & . & . \\ . & . & . & . & . & . & . & . & . & . & . \\ . & . & . & . & . & . & . & C_{n-2} & B_{n-2} & A_{n-2} & 0 & \\ . & . & . & . & . & . & . & 0 & C_{n-1} & B_{n-1} & A_{n-1} \\ . & . & . & . & . & . & . & . & 0 & 0 & 1 \\ \end{bmatrix} \quad \begin{bmatrix} y_0 \\ y_1 \\ y_2 \\ . \\ . \\ . \\ . \\ y_{n-1} \\ y_n \end{bmatrix} \qquad &= \qquad \begin{bmatrix} 1 \\ 4x_1^2 \delta x^2 \\ 4x_2^2 \delta x^2 \\ . \\ . \\ . \\ 4x_{n-1}^2 \delta x^2 \\ 6 \end{bmatrix} \\ \text{ } \\ A \qquad\qquad\qquad\qquad\qquad \boldsymbol{x} \qquad\quad &= \quad\qquad \boldsymbol{b} \end{align} $$




2.3.1 Code to solve numerically¶

Write functions to enable re-use and speed up iterations for multiple values of $n$.

In [110]:
def Aterm( partition_size: float) -> float:
    return 1 + (partition_size/2) * 3
    
def Bterm( partition_size: float) -> float:
    return -2 + (partition_size**2) * 2    

def Cterm( partition_size: float) -> float:
    return 1 - (partition_size/2) * 3

def fx( x_val: float) -> float:
    return 4 * (x_val ** 2)

def solve_diff( n: int, lower: int, upper: int, y0: float, y_inf: float,
                Afunc, Bfunc, Cfunc, funcf, debug: bool=False) -> np.array:
    dx = (upper - lower) /  (n+1)
    arr_x = np.linspace( lower, upper, n+1)
    arr_f = np.hstack( (y0, (dx**2) * np.vectorize(funcf)(arr_x[1:n]), y_inf))
    arr_k = (dx**2) * np.vectorize(funcf)(arr_x[1:n])
    arr_A = np.zeros( (n+1, n+1))
    # assign boundary conditions
    arr_A[0,0] = 1
    arr_A[n,n] = 1
    # assign A, B, C terms
    arr_terms = np.array([Cfunc(dx), Bfunc(dx), Afunc(dx)]) 
                # np.array([Afunc(dx), Bfunc(dx), Cfunc(dx)])
    for k in range(1,n):
        arr_A[k, k-1:k+2] = arr_terms
    
    if debug: 
        print(f"{dx=}\t{dx**2=}")
        print(f"{arr_x=}\n")
        print(f"{arr_f=}\n")
        print(f"{arr_k=}\n")
        print(f"{arr_A.shape=}\n")
        print(f"\narr_A=\n{arr_A}\n")
    arr_out = np.hstack(( arr_x.reshape(-1,1), 
                          np.linalg.solve( arr_A, arr_f).reshape(-1,1) ))
    return arr_out

def to_df( arr: np.array) -> pd.DataFrame:
    return pd.DataFrame( arr, columns=["x", "y"])

def style_table(df: pd.DataFrame, columns: int=3, spacing: int=20):
    space = "\xa0" * spacing
    k, m = divmod(df.shape[0], columns)
    output = ""
    for j in range(columns):
        df_table = df[j*k+min(j, m):(j+1)*k+min(j+1, m)]
        df_styler = df_table.style.set_table_attributes("style='display:inline'")
        output += df_styler._repr_html_() + space
    return output

For $n=10$,
In [111]:
df10 = to_df(solve_diff( 10, 1, 2, 1., 6., Aterm, Bterm, Cterm, fx))
display_html(style_table( df10, 1) , raw=True)
  x y
0 1.000000 1.000000
1 1.100000 2.299256
2 1.200000 3.288446
3 1.300000 4.034290
4 1.400000 4.591614
5 1.500000 5.005412
6 1.600000 5.312546
7 1.700000 5.543168
8 1.800000 5.721885
9 1.900000 5.868738
10 2.000000 6.000000
                    

For $n=50$,
In [112]:
df50 = to_df( solve_diff( 50, 1, 2, 1, 6, Aterm, Bterm, Cterm, fx))
display_html(style_table( df50) , raw=True)
  x y
0 1.000000 1.000000
1 1.020000 1.306593
2 1.040000 1.596244
3 1.060000 1.869767
4 1.080000 2.127942
5 1.100000 2.371518
6 1.120000 2.601211
7 1.140000 2.817710
8 1.160000 3.021674
9 1.180000 3.213736
10 1.200000 3.394503
11 1.220000 3.564556
12 1.240000 3.724452
13 1.260000 3.874727
14 1.280000 4.015892
15 1.300000 4.148439
16 1.320000 4.272837
                    
  x y
17 1.340000 4.389538
18 1.360000 4.498974
19 1.380000 4.601560
20 1.400000 4.697691
21 1.420000 4.787748
22 1.440000 4.872095
23 1.460000 4.951081
24 1.480000 5.025039
25 1.500000 5.094290
26 1.520000 5.159140
27 1.540000 5.219882
28 1.560000 5.276797
29 1.580000 5.330154
30 1.600000 5.380209
31 1.620000 5.427211
32 1.640000 5.471393
33 1.660000 5.512981
                    
  x y
34 1.680000 5.552192
35 1.700000 5.589231
36 1.720000 5.624296
37 1.740000 5.657576
38 1.760000 5.689252
39 1.780000 5.719495
40 1.800000 5.748471
41 1.820000 5.776338
42 1.840000 5.803246
43 1.860000 5.829340
44 1.880000 5.854756
45 1.900000 5.879627
46 1.920000 5.904078
47 1.940000 5.928229
48 1.960000 5.952195
49 1.980000 5.976083
50 2.000000 6.000000
                    

For $n=100$,

In [113]:
df100 = to_df( solve_diff( 100, 1, 2, 1, 6, Aterm, Bterm, Cterm, fx))
display_html( style_table( df100, 4, 5), raw=True)
  x y
0 1.000000 1.000000
1 1.010000 1.156879
2 1.020000 1.309337
3 1.030000 1.457482
4 1.040000 1.601420
5 1.050000 1.741253
6 1.060000 1.877083
7 1.070000 2.009009
8 1.080000 2.137128
9 1.090000 2.261535
10 1.100000 2.382323
11 1.110000 2.499583
12 1.120000 2.613404
13 1.130000 2.723873
14 1.140000 2.831077
15 1.150000 2.935098
16 1.160000 3.036018
17 1.170000 3.133918
18 1.180000 3.228876
19 1.190000 3.320969
20 1.200000 3.410273
21 1.210000 3.496860
22 1.220000 3.580803
23 1.230000 3.662172
24 1.240000 3.741037
25 1.250000 3.817465
     
  x y
26 1.260000 3.891523
27 1.270000 3.963274
28 1.280000 4.032783
29 1.290000 4.100111
30 1.300000 4.165320
31 1.310000 4.228469
32 1.320000 4.289615
33 1.330000 4.348816
34 1.340000 4.406128
35 1.350000 4.461605
36 1.360000 4.515301
37 1.370000 4.567267
38 1.380000 4.617555
39 1.390000 4.666215
40 1.400000 4.713296
41 1.410000 4.758846
42 1.420000 4.802911
43 1.430000 4.845538
44 1.440000 4.886771
45 1.450000 4.926655
46 1.460000 4.965232
47 1.470000 5.002544
48 1.480000 5.038632
49 1.490000 5.073537
50 1.500000 5.107299
     
  x y
51 1.510000 5.139954
52 1.520000 5.171542
53 1.530000 5.202099
54 1.540000 5.231662
55 1.550000 5.260264
56 1.560000 5.287942
57 1.570000 5.314728
58 1.580000 5.340656
59 1.590000 5.365757
60 1.600000 5.390065
61 1.610000 5.413608
62 1.620000 5.436418
63 1.630000 5.458525
64 1.640000 5.479956
65 1.650000 5.500741
66 1.660000 5.520906
67 1.670000 5.540480
68 1.680000 5.559488
69 1.690000 5.577956
70 1.700000 5.595909
71 1.710000 5.613372
72 1.720000 5.630370
73 1.730000 5.646926
74 1.740000 5.663062
75 1.750000 5.678802
     
  x y
76 1.760000 5.694167
77 1.770000 5.709180
78 1.780000 5.723861
79 1.790000 5.738230
80 1.800000 5.752308
81 1.810000 5.766115
82 1.820000 5.779670
83 1.830000 5.792991
84 1.840000 5.806097
85 1.850000 5.819006
86 1.860000 5.831735
87 1.870000 5.844302
88 1.880000 5.856723
89 1.890000 5.869015
90 1.900000 5.881193
91 1.910000 5.893273
92 1.920000 5.905271
93 1.930000 5.917202
94 1.940000 5.929079
95 1.950000 5.940917
96 1.960000 5.952730
97 1.970000 5.964532
98 1.980000 5.976336
99 1.990000 5.988154
100 2.000000 6.000000
     



2.3.2 Closed form general solution of ODE¶

Obtain the closed form solution of the ODE using SymPy.

In [114]:
y = sp.symbols("y", cls=sp.Function)
x = sp.symbols("x")
eqs = y(x).diff(x).diff(x) + 3*y(x).diff(x) + 2*y(x) - 4*x**2

# boundary condition
bound_cond = {y(1):1, y(2):6}    
general_soln = sp.dsolve(eqs, ics=bound_cond)
general_soln
Out[114]:
$\displaystyle y{\left(x \right)} = 2 x^{2} - 6 x + 7 + \frac{\left(2 e + 3 e^{3}\right) e^{- x}}{-1 + e} + \frac{\left(- 3 e^{4} - 2 e^{3}\right) e^{- 2 x}}{-1 + e}$



Generate dataframe for plotting.

In [115]:
xval = np.linspace(1,2,600)
y_func = sp.lambdify( x, general_soln.rhs)
df_closedform = pd.DataFrame( np.hstack((xval.reshape(-1,1), y_func(xval).reshape(-1,1) )),
                              columns=['x','y'])
In [116]:
df_closedform
Out[116]:
x y
0 1.000000 1.000000
1 1.001669 1.026755
2 1.003339 1.053382
3 1.005008 1.079881
4 1.006678 1.106253
... ... ...
595 1.993322 5.992169
596 1.994992 5.994125
597 1.996661 5.996082
598 1.998331 5.998041
599 2.000000 6.000000

600 rows × 2 columns


Generate numerical solution for varying $n$ and plot. (Values of $n$ chosen for visual aesthetics)

In [117]:
ls_df = []
n_vals = [i for i in range(2,10,2)]
for k in n_vals:
    ls_df.append( to_df( solve_diff( k, 1, 2, 1, 6, Aterm, Bterm, Cterm, fx)))
In [118]:
fig12 = px.scatter( title="<b>Figure 12 - Varying size of <i>n</i></b>", width=850, height=600
                ).update_layout(yaxis_range=[0.9,6.1],
                                xaxis_range=[0.95,2.05]) 

for i in range(len(ls_df)):
    fig12.add_scatter(mode='lines+markers', name=f"n={n_vals[i]:,}",
                    x = ls_df[i].x,
                    y = ls_df[i].y,
                    hovertemplate = 'x: <b>%{x:.6f}</b><br>y: <b>%{y:,.4f}</b>',
                    line=dict( width=1.5), marker=dict(symbol='triangle-down', size=10)
                   )

fig12.add_scatter(mode='lines', name=f"Closed Form",
                x = df_closedform.x,
                y = df_closedform.y,
                hovertemplate = 'x: <b>%{x:.6f}</b><br>y: <b>%{y:,.4f}</b>',
                line=dict( width=2.5, color='black', dash='dash')
               )
# Show spikes
fig12.update_xaxes(showspikes=True, title='x', dtick=0.1)
fig12.update_yaxes(showspikes=True, title='y')
fig12.show()

As $n$ increases, the numerical solution approaches the exact solution. Essentially, $n$ is the number of straight lines to approximate a curve, the larger the value, the finer and better able it is match the exact solution.





2.3.3 Plot the error percentage¶

Let error be the difference between numerical and general solution.

$$ \varepsilon = \hat{y} - y \\ $$
In [119]:
ls_df = []
n1_vals = [2, 5, 10, 25, 50, 100, 200, 500, 1000, 2000, 5000] 
for k in n1_vals:
    closed_form = y_func( np.linspace(1,2,k+1))
    df_num_soln = to_df( solve_diff( k, 1, 2, 1, 6, Aterm, Bterm, Cterm, fx))
    df_num_soln['exact'] = closed_form
    df_num_soln['error'] = df_num_soln.y - df_num_soln.exact
    df_num_soln['error_percentage'] = 100 * (df_num_soln.error / closed_form[1])
    ls_df.append( df_num_soln )
In [120]:
ls_df[0].tail()
Out[120]:
x y exact error error_percentage
0 1.0 1.00000 1.000000 0.000000 0.000000
1 1.5 4.78125 5.120806 -0.339556 -6.630906
2 2.0 6.00000 6.000000 0.000000 0.000000
In [121]:
fig13 = px.scatter( title="<b>Figure 13 - Error percentage for varying size of <i>n</i></b>",
                   width=750, height=800)

for k in range(len(ls_df)):
    fig13.add_scatter(mode='lines', name=f"n={n1_vals[k]:,}",
                    x = ls_df[k].x,
                    y = ls_df[k].error_percentage,
                    hovertemplate = 'x: <b>%{x:.6f}</b><br>Error (%): <b>%{y:.5f}</b>',
                    line=dict( width=1.5))

# Show spikes
fig13.update_xaxes(showspikes=True, title='x', dtick=0.1)
fig13.update_yaxes(showspikes=True, title='Error (%)', tickformat='0.1f')
fig13.show()

Based on the chart above, at around 200 partitions of $x$ the maximum error is less than 1%. And at 2000 partitions, less than 0.1% error. Note that $x=1$ and $x=2$ are the boundary conditions and are fixed.
Negative indicates that the numerical estimate is underestimating the solution.











Question 3¶

Solve the following using Monte Carlo integration

$$ \int_{1}^{3} x^2 \,dx \tag{Integral 1} \\[10pt] $$$$ \int_{0}^{\infty} e^{-x^2} \,dx \tag{Integral 2} \\[10pt] $$$$ \frac {1}{\sqrt{2 \pi}} \int_{-\infty}^{\infty} x^4 e^{-x^2/2} \,dx \tag{Integral 3} \\[10pt] $$

3.1 Description¶

Monte Carlo method is centered on evaluating definite integrals as expectations.

For example, to estimate $\theta$ given that $$ \theta = \mathbb{E}[f(U)] \quad \text{where} \quad U \sim U(0,1) $$

The expectation can be expressed as an integral of the form $$ \mathbb{E}[f(U)] = \int_{0}^{1} f(u) p(u) \,du $$

where $p(u)$ is the probability density function (PDF) for the uniform distribution $U(,0,1)$ given by

$$ p(u) = \begin{cases} 1 & 0 \le u \leq1 \\ 0 & \mbox{otherwise} \end{cases} $$

hence $$ \mathbb{E}[f(U)] = \int_{0}^{1} f(u) \,du $$


There are two fundamental cases depending on the limits of integration.

Case 1:   function of $f$   $\text{s.t.} \quad f:[a,b] \rightarrow \mathbb{R} \quad \text{ where } \quad -\infty \lt a \lt b \lt \infty$

For a given integral, $$ I = \int_{a}^{b} f(x)\,dx \quad \\[10pt] $$

transform the domain $[a,b]$ to $[0,1]$ by applying the subsititution,

$$ y = \frac{x-a}{b-a} \quad \Rightarrow \quad x = y \text{ } (b-a) + a $$

which gives

$$ dx = dy \text{ } (b-a) $$

So now,

$$ \begin{align} I \text{ } &= \text{ } (b-a) \int_{0}^{1} f \left( \text{ } y(b-a) + a \text{ } \right) \,dy \\[8pt] &= \text{ } (b-a) \text{ }\text{ } \mathbb{E} \left[ \text{ } f( \text{ } U (b-a) + a \text{ } ) \text{ } \right] \quad \text{where} \quad U \sim U(0,1) \end{align} $$



Case 2:   function of $g$   $\text{s.t.} \quad g:[a,\infty) \rightarrow \mathbb{R} \text{ where } \quad -\infty \lt a \lt b \lt \infty$

For a given integral $I$, provided $I \lt \infty$ $$ I = \int_{a}^{b} g(x)\,dx \\[10pt] $$

transform the domain $[a,b]$ to $[0,1]$ by applying the subsititution,

$$ y = \frac{1}{1+x} \quad \Rightarrow \quad x = \frac{1}{y} - 1 $$

which gives

$$ dx = \frac{dy} {y^2} $$

So now,

$$ \begin{align} I &= \int_{0}^{1} \frac {g \left( \text{ } \frac{1}{y} - 1 \text{ } \right)} {y^2} \,dy \\[10pt] &= \mathbb{E} \left[ \frac{ g \left( -1 + \frac{1}{U} \right) } { U^2} \right] \quad \text{where} \quad U \sim U(0,1) \end{align} $$



Dr Ahmad, R. (2023). Monte Carlo Tutorial - "Monte Carlo" [PowerPoint slides].
CQF Program January 2023

In all proceeding cases, define error to be the difference between estimate and exact. $$ \varepsilon = \hat{I_1} - I_1 $$


3.2 Integral 1¶

3.2.1 Exact solution¶

Given the integral,

$$ I_1 = \int_{1}^{3} x^2 \,dx $$

the exact solution is,

$$ \begin{align} I_1 &= \left[\frac{x^3}{3}\right]_1^3 \;=\; \left[\frac{3^3}{3}\right] - \left[\frac{1^3}{3}\right] \\[7pt] &= \left[\frac{27}{3}\right] - \left[\frac{1}{3}\right] \\[7pt] &= \frac{26}{3} \end{align} $$

3.2.2 Monte Carlo Integration¶

Let,
$$ \begin{align} I_1 &= \int_{1}^{3} f(x) \,dx \quad \text{where} \quad f(x) = x^2 \\[7pt] \end{align} $$

applying Case 1 to transform from $[1,3]$ to domain $[0,1]$ by the following substitution,

$$ y = \frac{x-1}{3-1} = \frac{x-1}{2} \quad \Rightarrow \quad x = 2y + 1 $$

when
 $ x = 1$,  $ y = 0$
 $ x = 3$,  $ y = 1$


differentiating $$ dx = 2 \text{ } dy $$

This gives,

$$ \begin{align} I_1 \text{ } &= \text{ } \int_{0}^{1} f( \text{ } 2 \text{ } y + 1 \text{ } ) \, \text{ } 2 dy \;=\; \text{ } 2 \text{ } \int_{0}^{1} f( \text{ } 2 \text{ } y + 1 \text{ } ) \,dy \\[5pt] &= \text{ } 2 \text{ }\text{ } \mathbb{E} \left[ \text{ } f( \text{ } 2U+1 \text{ } ) \text{ } \right] \quad \text{where} \quad U \sim U(0,1) \end{align} $$

Next,

  1. Begin by generating a sequence $U_i \sim U(0,1)$
  2. Compute $Y_i = (2U_i+1)^2 $ for $ i = 1,...,n$
  3. Calculate the estimate $ \hat{I_1}$
$$ \begin{align} \hat{I_1} &= 2 \text{ } \left( \frac{1}{n} \sum_{1}^{n} Y_i \right) \\ &= 2 \text{ } \left( \frac{1}{n} \sum_{1}^{n} ( \text{ } 2U_i+1 \text{ })^2 \right) \end{align} $$
In [122]:
exact_I1 = 26 / 3
iters_arr1 = 10**7

def solve_I1(exact: float, iter_arr1: int) -> pd.DataFrame:
    np.random.seed(20230414)
    arr1 = np.random.uniform(0, 1, iter_arr1)
    arr1 = np.vstack((arr1, (arr1*2+1)**2))                 # 1 - term
    arr1 = np.vstack((arr1, (np.cumsum(arr1[1,:]) / np.arange(1,arr1.shape[1]+1)) ))    # 2 - average
    # extract values at interval - save RAM and swap
    arr1_interval = np.append(1, np.arange(1_000,arr1.shape[1], 1_000))
    arr1 = arr1[:,arr1_interval]
    arr1 = np.vstack((arr1, arr1[2,:] * 2))                 # 3 - average*2
    arr1 = np.vstack((arr1, arr1[3,:] - exact))             # 4 - error            
    arr1 = np.vstack((arr1, (arr1[4,:] / exact) * 100 ))    # 5 - error percentage  
    arr1 = np.transpose(np.vstack((arr1, arr1_interval)))
    df = pd.DataFrame( arr1, 
                      columns=['rand','term','average', 'estimate', 'error', 'error_perc','iter'])
    return df
In [123]:
df1 = solve_I1(exact_I1, iters_arr1)
# df1.shape
print(f"The closed form solution of Integal 1 has the value {exact_I1:0.6f}\n")
The closed form solution of Integal 1 has the value 8.666667

In [124]:
df1.iloc[:,3:].tail()
Out[124]:
estimate error error_perc iter
9995 8.666773 0.000107 0.001230 9995000.0
9996 8.666760 0.000093 0.001078 9996000.0
9997 8.666764 0.000097 0.001125 9997000.0
9998 8.666788 0.000121 0.001401 9998000.0
9999 8.666793 0.000127 0.001460 9999000.0

Table above shows part of the calculated values for Integral 1. Compared to the exact solution, the estimate appears to have converged.


In [125]:
arr1_fig1 = px.scatter( df1, x='iter', y='estimate',
                 labels={'estimate': 'Estimate', 'iter': 'Iterations'},
                 title="<b>Figure 14 - Estimate vs Iteration for Integral 1</b><br><sup>(exact solution in red)</sup>",
                 width=900, height=500,
                ).update_traces(mode='lines', line = dict(color='blue'),
                                hovertemplate='Error: %{y}<br>Iteration: %{x}'
                ).update_layout(yaxis_range=[8.64,8.68])
arr1_fig1.add_hline(y=exact_I1, line_width=1.5, line_dash="solid", line_color="red", opacity=0.7, name="exact")
# Show spikes
arr1_fig1.update_xaxes(showspikes=True)
arr1_fig1.update_yaxes(showspikes=True, tickformat='.2f')
arr1_fig1.show()

Chart above shows the convergence of the Monte Carlo integration for Integral 1. As can be observed, the numerical integration converges relatively fast to some small error after some initial fluctuations. The estimate never seems to truly converge and appears to oscillate around the closed form value. The next chart will show the error percentage to have better sense of performance.


In [126]:
arr1_fig3 = px.scatter( df1, x='iter', y='error_perc',
                 labels={'error_perc': 'Error (%)', 'iter': 'Iterations'},
                 title="<b>Figure 15 - Iteration vs Error (%) for Integral 1</b>",
                 width=900, height=600,
                ).update_traces(mode='lines', line = dict(color='red'),
                                hovertemplate='Error (%): %{y:.6%}<br>Iteration: %{x}'                                
                ).update_layout(yaxis_range=[-0.10,0.06])
arr1_fig3.add_hline(y=0, line_width=1, line_dash="solid", line_color="black", opacity=0.3)
# Show spikes
arr1_fig3.update_xaxes(showspikes=True)
arr1_fig3.update_yaxes(showspikes=True, tickformat='.2%', dtick = 0.01)
arr1_fig3.update_layout(showlegend=True)
arr1_fig3.show()

The plot above shows the error as percentage against iterations for the Monte Carlo integration. After approximately 7 million iterations, the error is within 1% range. Depending on use case, this may or may not be sufficient accuracy.



3.3 Integral 2¶

3.3.1 Exact solution¶

Given, $$ \begin{align} I_2 &= \int_{0}^{\infty} e^{-x^2} \,dx \tag{Integral 2} \\[10pt] &= \frac{1}{2} \text{ } \underbrace{ \int_{-\infty}^{\infty} e^{-x^2} \,dx }_{\text{Gaussian integral}} \end{align}\\ $$ since it is known that the Gaussian interval is symmetric about $x=0$.

Notice the Gaussian integral which the answer is known

$$ {\int_{-\infty}^{\infty} e^{-x^2} \,dx} = \sqrt{\pi} $$

Therefore $$ I_2 = \frac{\sqrt{\pi}} {2} $$

3.3.2 Monte Carlo Integration¶

Let,

$$ I_2 = \int_{0}^{\infty} g(x) \,dx \quad \text{where} \quad g(x) = e^{-x^2} \tag{Integral 2} \\[10pt] $$

applying Case 2 to transform from $[0,\infty]$ to domain $[0,1]$ by the following substitution,

$$ y = \frac{1}{1 + x} \quad \Rightarrow \quad x = \frac{1}{y} - 1 $$

when
 $ x = 0$,    $ y = 1$
 $ x = \infty$,  $ y = 0$


and differentiating gives

$$ dx = - \text{ } \frac{dy} {y^2} $$

This gives,

$$ \begin{align} I_1 \text{ } &= - \int_{1}^{0} g( \frac{1}{y} - 1 ) \text{ } \, \frac{ dy }{y^2} \\[7pt] &= \int_{0}^{1} \frac{ g( \frac{1}{y} - 1 )}{y^2} \text{ } \, dy \\[7pt] &= \text{ } \mathbb{E} \left[ \frac{ g( \frac{1}{U}-1) }{U^2} \right] \quad \text{where} \quad U \sim U(0,1) \end{align} $$

Next,

  1. Begin by generating a sequence $U_i \sim U(0,1)$
  2. Compute $ \quad Y_i = \frac{ e^{ - \text{ }( \frac{1}{U_i} - 1)^2}}{U_I^2} \quad \text{for} \quad i = 1,...,n $
  3. Calculate the estimate $ \hat{I_2}$
$$ \begin{align} \hat{I_2} &= \text{ } \frac{1}{n} \sum_{1}^{n} Y_i \\ &= \text{ } \frac{1}{n} \sum_{1}^{n} \frac{ e^{ - \text{ }( \frac{1}{U_i} - 1)^2}} {U_I^2} \end{align} $$
In [127]:
exact_i2 = np.sqrt(np.pi) / 2
iter_arr2 = 10**7

def solve_I2( exact2: float, iters: int) -> pd.DataFrame:
    np.random.seed(20230414)
    arr2 = np.random.uniform(0, 1, iters)
    arr2 = np.vstack( (arr2, arr2**2))                              # 1 - denominator
    arr2 = np.vstack( (arr2, np.exp( -((1/arr2[0,:]) - 1)**2) ))    # 2 - numerator
    arr2 = np.vstack( (arr2, arr2[2,:] / arr2[1,:]))                # 3 - numerator / denominator
    arr2 = np.vstack( (arr2, np.cumsum(arr2[3,:]) / np.arange(1,arr2.shape[1]+1) ))    # 4- estimate (averaged) 
    # extract values at interval
    arr2_interval = np.append(1, np.arange(2_000,arr2.shape[1], 2_000))
    arr2 = arr2[:,arr2_interval]
    arr2 = np.vstack( (arr2, arr2[4,:] - exact2))                   # 5 - error
    arr2 = np.vstack( (arr2, (arr2[5,:] / exact2) * 100))           # 6 - error percentage
    arr2_plot = np.transpose( np.vstack( (arr2, arr2_interval )))
    df = pd.DataFrame( arr2_plot,
                       columns=['randuniform', 'denom', 'numer', 'term', 'estimate','error',
                                'error_percentage','iter'])
    return df
In [128]:
df2 = solve_I2(exact_i2, iter_arr2)
print(f"The closed form solution of Integal 2 has the value {exact_i2:0.6f}\n")
The closed form solution of Integal 2 has the value 0.886227

In [129]:
df2.iloc[:,4:].tail()
Out[129]:
estimate error error_percentage iter
4995 0.886315 0.000088 0.009919 9990000.0
4996 0.886318 0.000091 0.010221 9992000.0
4997 0.886318 0.000091 0.010234 9994000.0
4998 0.886317 0.000090 0.010157 9996000.0
4999 0.886318 0.000091 0.010297 9998000.0

Table above shows part of the calculated values for Integral 2. Compared to the exact value, the estimate appears to have converged.


In [130]:
arr2_fig1 = px.scatter( df2, x='iter', y='estimate',
                 labels={'estimate': 'Estimate', 'iter': 'Iterations'},
                 title="<b>Figure 16 - Estimate vs Iteration for Integral 2</b><br><sup>(exact solution in red)</sup>",
                 width=900, height=500,
            ).update_traces(mode='lines', hovertemplate='Estimate: %{y}<br>Iteration: %{x}'
            ).update_layout(yaxis_range=[0.882,0.888])
arr2_fig1.add_hline(y=exact_i2, line_width=1, line_dash="solid", line_color="red", opacity=0.9)
# Show spikes
arr2_fig1.update_xaxes(showspikes=True)
arr2_fig1.update_yaxes(showspikes=True, tickformat=".3f")
arr2_fig1.show()

Figure above show the estimate performance for increasing iterations. The estimate appears to converge beyond 6 million iterations with some minor oscillations around the exact solution.


In [131]:
arr2_fig3 = px.scatter( df2, x='iter', y='error_percentage',
                 labels={'error_percentage': 'Error (%)', 'iter': 'Iterations'},
                 title="<b>Figure 17 - Error (%) vs Iteration for Integral 2</b>",
                 width=900, height=500,
                ).update_traces(mode='lines', line = dict(color='red'),
                                hovertemplate='Error (%): %{y}<br>Iteration: %{x}'
                ).update_layout(yaxis_range=[-0.1,0.1])
arr2_fig3.add_hline(y=0, line_width=1, line_dash="solid", line_color="black", opacity=0.3)
arr2_fig3.add_hline(y=0.015, line_width=1.5, line_dash="dash", line_color="green", opacity=1)
arr2_fig3.add_hline(y=-0.015, line_width=1.5, line_dash="dash", line_color="green", opacity=1)
# Show spikes
arr2_fig3.update_xaxes(showspikes=True)
arr2_fig3.update_yaxes(showspikes=True, tickformat='.2%', dtick=0.02) 
arr2_fig3.show()

The plot above shows the error between the estimate and the exact solution. It appears to be fluctuating up until 5 million iterations before starting to converge and even after that, it does not seem to stay within $\pm$1.5% (green dashed line) and oscillates. Note that for the Gaussian definite integral, a solution is only possible using multivariable calculus, a numerical answer with approximately 1.5% error after 10 million iterations is perhaps passable depending on the situation. Unfortunately, there are no guarantees that more iterations will finally converge. It will be left for further investigation.



3.4 Integral 3¶

3.4.1 Exact solution¶

Given the integral, $$ \begin{align} I_3 &= \frac {1}{\sqrt{2 \pi}} \int_{-\infty}^{\infty} x^4 e^{-x^2/2} \,dx \\[5pt] \end{align} $$

Apply integration by parts repeatedly onto the integral,

$$ \begin{align} \int_{-\infty}^{\infty} x^4 e^{-x^2/2} \,dx \quad &= \quad \int_{-\infty}^{\infty} x^3 \text{ } \mbox{ } x \text{ } e^{-x^2/2} \,dx \\[10pt] &= \quad \underbrace{ \left[ - e^{-x^2/2} \text{ } x^3 \right]_{-\infty}^{\infty} }_{\lim_{ x\to\infty} \rightarrow 0 } + 3 \int_{-\infty}^{\infty} x^2 e^{-x^2/2} \,dx \\[10pt] &= 3 \text{ } \left( \text{ } \underbrace{ \left[ -e^{-x^2/2} \text{ } x \right]_{-\infty}^{\infty} }_{\lim_{ x\to\infty} \rightarrow 0 } + \underbrace{ \int_{-\infty}^{\infty} e^{-x^2/2} \,dx }_{{\textstyle\unicode{x24B6}}} \text{ } \right) \\[10pt] &= 3 \sqrt{2 \pi} \\[5pt] \end{align} $$

For the integral $\textstyle\unicode{x24B6}$,

$$ I = \int_{-\infty}^{\infty} e^{-x^2/2} \,dx = \int_{-\infty}^{\infty} e^{-y^2/2} \,dy \\[15pt] $$

Then,

$$ \begin{align} I^2 &= \int_{-\infty}^\infty \int_{-\infty}^\infty e^{ -(x^2+y^2)/2 } \,dx \,dy \\[5pt] &= \int_{0}^{2\pi} \int_{0}^{\infty} e^{-r^2/2} \text{ } r \,dr \,d\theta \\[5pt] &= 1 \text{ } \int_0^{2\pi} \,d\theta \\[5pt] &= 2\pi \\[10pt] \Rightarrow \quad I &= \sqrt{ 2 \pi} \\[5pt] \end{align}\\ $$

The exact solution is $$ \begin{align}\\ I_3 &= \frac {1}{\sqrt{2 \pi}} \text{ } \left( \text{ } 3 \sqrt{2 \pi} \text{ } \right) \\[5pt] &= 3 \\[5pt] \end{align}\\ $$

3.4.2 Monte Carlo Integration¶

Let,

$$ \begin{align} I_3 &= \frac {1}{\sqrt{2 \pi}} \int_{-\infty}^{\infty} g(x) \,dx \\[10pt] &= \frac {1}{\sqrt{2 \pi}} \text{ } \left[ \text{ } \underbrace{ \int_{-\infty}^{0} g(x) \,dx }_{{\textstyle\unicode{x2460}}}\text{ } + \text{ } \underbrace{ \int_{0}^{\infty} g(x) \,dx }_{{\textstyle\unicode{x2461}}} \text{ } \right] \\[12pt] & \text{where} \quad g(x) = x^4 e^{-x^2/2} \end{align} $$

applying Case 2 to transform from $[0,\infty]$ to domain $[0,1]$ by the following substitution,


$$ y = \frac{1}{1 + x} \quad \Rightarrow \quad x = \frac{1}{y} - 1 $$

for integral $ \textstyle\unicode{x2460} $ applying the transformation yields the new domain
 $ x = -\infty $,    $ y = 0$
 $ x = 0 $,    $ y = 1$

and for integral $ \textstyle\unicode{x2461} $
 $ x = 0$,    $ y = 1$
 $ x = \infty$,  $ y = 0$


and differentiating gives

$$ dx = - \text{ } \frac{dy} {y^2} $$

This gives,

$$ \begin{align} I_3 \text{ } &= \frac {1}{\sqrt{2 \pi}} \text{ } \left( \text{ } \int_{0}^{1} g( \frac{1}{y} - 1 ) \text{ } \, \frac{dy}{y^2} \quad - \quad \int_{1}^{0} g( \frac{1}{y} - 1 ) \text{ } \, \frac{dy}{y^2} \text{ } \right) \\[10pt] &= \frac {1}{\sqrt{2 \pi}} \text{ } \left( \text{ } \int_{0}^{1} \frac{ g( \frac{1}{y} - 1) }{y^2} \text{ } \, dy \quad + \quad \int_{0}^{1} \frac{ g( \frac{1}{y} - 1) }{y^2} \text{ } \, dy \text{ } \right) \\[10pt] &= \frac {2}{\sqrt{2 \pi}} \text{ } \left( \text{ } \int_{0}^{1} \frac{ g( \frac{1}{y} - 1) }{y^2} \text{ } \, dy \text{ } \right) \\[10pt] &= \frac {2}{\sqrt{2 \pi}} \text{ } \mathbb{E} \left[ \frac{ g( \frac{1}{U}-1) }{U^2} \right] \quad \text{where} \quad U \sim U(0,1) \end{align} $$

Next,

  1. Begin by generating a sequence $U_i \sim U(0,1)$
  2. Compute $ \quad Y_i = ( \frac{1}{U_i} - 1)^4 \text{ } e^{ \frac{ - \text{ }( \frac{1}{U_i} - 1)^2 }{2} } \text{ } U_i^{-2} \quad \text{for} \text{ } i = 1,...,n $
  3. Calculate the estimate $ \hat{I_3}$
$$ \begin{align} \hat{I_3} &= \text{ } \frac {2}{\sqrt{2 \pi}} \text{ } \left( \frac{1}{n} \sum_{1}^{n} Y_i \right) \\[7pt] &= \text{ } \frac {2}{\sqrt{2 \pi}} \text{ } \left( ( \frac{1}{U_i} - 1)^4 \text{ } e^{ \frac{ - \text{ }( \frac{1}{U_i} - 1)^2 }{2} } \text{ } U_i^{-2} \right) \\[7pt] \end{align} $$
In [132]:
exact_i3 = 3
iter_arr3 = 10**8

def solve_I3( exact: float, iter: int) -> pd.DataFrame:
    np.random.seed(20230414)
    arr3 = np.random.uniform(0, 1, iter)
    arr3 = np.vstack(( arr3, ((1/arr3) - 1) ))                              # 1 - input to function
    arr3 = np.vstack(( arr3, np.square(arr3[0,:]) ))                        # 2 - U^2
    arr3 = np.vstack(( arr3, arr3[1,:]**4 ))                                # 3 - x^4 (first term)
    arr3 = np.vstack(( arr3, np.exp(-np.square(arr3[1,:]) / 2) ))           # 4 - exp (second term)
    arr3 = np.vstack(( arr3, arr3[3,:] * arr3[4,:] * (1/ arr3[2,:]) ))      # 5 - all 
#     arr3 = np.vstack(( arr3, np.cumsum(arr3[5,:]) / np.arange(1,arr3.shape[1]+1) ))    # 6 - averaged 
#     arr3 = np.vstack(( arr3, (2/np.sqrt(np.pi*2)) * arr3[6,:] ))            # 7 - estimate
    arr3 = np.vstack(( arr3, (2/np.sqrt(np.pi*2)) * arr3[5,:] ))            # 7 - estimate
    arr3 = np.vstack(( arr3, np.cumsum(arr3[6,:]) / np.arange(1,arr3.shape[1]+1) ))    # 6 - averaged 
    
    # extract values at interval
    arr3_interval = np.append(1, np.arange(10_000, arr3.shape[1], 10_000))
    arr3 = arr3[:, arr3_interval]
    arr3 = np.vstack(( arr3, arr3[7,:] - exact ))                          # 8 - error
    arr3 = np.vstack(( arr3, (arr3[8,:] / exact) * 100 ))                  # 9 - error percentage
    arr3_plot = np.transpose( np.vstack( (arr3, arr3_interval )))
    df = pd.DataFrame( arr3_plot,
                       columns=['randuniform', 'input_arg', 'u_square', 'first_term', 'second_term',
                                'all_terms','averaged', 'estimate', 'error','error_percentage','iter'])
    return df
In [133]:
df3 = solve_I3( exact_i3, iter_arr3)
print(f"The closed form solution of Integal 2 has the value {exact_i3:0.6f}")
The closed form solution of Integal 2 has the value 3.000000
In [134]:
df3.iloc[:,7:].tail()
Out[134]:
estimate error error_percentage iter
9995 2.999926 -0.000074 -0.002465 99950000.0
9996 2.999933 -0.000067 -0.002225 99960000.0
9997 2.999932 -0.000068 -0.002278 99970000.0
9998 2.999933 -0.000067 -0.002231 99980000.0
9999 2.999930 -0.000070 -0.002323 99990000.0

Table above shows part of the calculated values for Integral 3. Compared to the exact value, the estimate appears to have converged.


In [135]:
arr3_fig1 = px.scatter( df3, x='iter', y='estimate',
                 labels={'estimate': 'Estimate', 'iter': 'Iterations'},
                 title="<b>Figure 18 - Estimate vs Iteration for Integral 3</b><br><sup>(exact solution in red)</sup>",
                 width=900, height=500,
            ).update_traces(mode='lines', hovertemplate='Estimate: %{y}<br>Iteration: %{x}'
            ).update_layout(yaxis_range=[2.990,3.003])
arr3_fig1.add_hline(y=exact_i3, line_width=1, line_dash="solid", line_color="red", opacity=0.5)
# Show spikes
arr3_fig1.update_xaxes(showspikes=True)
arr3_fig1.update_yaxes(showspikes=True, dtick=0.002, tickformat=".3f")
arr3_fig1.show()

Figure above show the estimate performance for increasing iterations. The estimate appears to converge beyond 6 million iterations with some minor oscillations around the exact solution.


In [136]:
arr3_fig3 = px.scatter( df3, x='iter', y='error_percentage',
                 labels={'error_percentage': 'Error (%)', 'iter': 'Iterations'},
                 title="<b>Figure 19 - Error (%) vs Iteration for Integral 3</b>",
                 width=900, height=500,
                ).update_traces(mode='lines', line = dict(color='red'),
                                hovertemplate='Error (%): %{y}<br>Iteration: %{x}'
                ).update_layout(yaxis_range=[-0.10,0.05])
arr3_fig3.add_hline(y=0, line_width=1, line_dash="solid", line_color="black", opacity=0.3)
arr3_fig3.add_hline(y=0.015, line_width=1, line_dash="dash", line_color="green", opacity=0.9)
arr3_fig3.add_hline(y=-0.015, line_width=1, line_dash="dash", line_color="green", opacity=0.9)
# Show spikes
arr3_fig3.update_xaxes(showspikes=True)
arr3_fig3.update_yaxes(showspikes=True, tickformat='.1%', dtick=0.02)
arr3_fig3.show()

The plot above shows the error between the estimate and the exact solution. Integral 3 required significantly higher iterations before the error was even within 1.5%. It appears to be fluctuating up until 50 million iterations before starting to converge to within $\pm$1.5% (green dashed line).



3.5 Conclusion¶

Monte Carlo integration is a very useful tool to approximate integrals and with the use of the numpy library, only takes 1-2 seconds to calculate depending on the iterations required and the accuracy desired. It is important to note that all 3 integrals seem to converge to within $\pm$1-2% of the exact solution. However, if this was used as part of a larger calcultions, the error accumulation has to be taken into account. Obviously, there is a speed and accuracy trade off. In the numerical estimations, no analysis has been done on stability of the solution. This should be part of the requirement when developing programs to frequently run as some integrals may never converge or contain asymptotic behavior over a certain range.

Further optimization can definitely be done to speed up the code including parallelization and reducing variable assignments.